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Abstract. We propose a novel approach to the problem of polynomial approximation of rational Bezier 
triangular patches with prescribed boundary control points. The method is very efficient thanks to 
using recursive properties of the bivariate dual Bernstein polynomials and applying a smart algorithm 
for evaluating a collection of two-dimensional integrals. Some illustrative examples are given. 
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1 Introduction and preliminaries 

Rational triangular Bezier surfaces are an important tool in computer graphics. However, 
they may be sometimes inconvenient in practical applications. The reason is that evaluation 
of integrals or derivatives of rational expressions is cumbersome. Also, it happens that a 
rational surface produced in one CAD system is to be imported into another system which 
can handle only polynomial surfaces. 

In order to solve the two problems above, different algorithms for approximating the ra¬ 
tional surface by polynomial surface are proposed [3,8,17,20-22]. The spectrum of methods 
contains hybrid algorithm [22], progressive iteration approximation [3,8], least squares ap¬ 
proximation and linear programming [17], and approximation by Bezier surfaces with control 
points obtained by successive degree elevation of the rational Bezier surface [20,21]. As a 
rule, no geometric constraints are imposed, which means a serious drawback: if we start 
with a patchwork of smoothly connected rational Bezier triangles and approximate each patch 
separately, we do not obtain a smooth composite surface. 

In this paper, we propose a method for solving the problem of the constrained least squares 
approximation of a rational triangular Bezier patch by a polynomial triangular Bezier patch; 
see Problem 2.1 below. The method is based on the idea of using constrained dual bivariate 
Bernstein polynomials. Using a fast recursive scheme of evaluation of Bezier form coefficients 
of the bivariate dual Bernstein polynomials, and applying a swift adaptive scheme of numerical 
computation of a collection of double integrals involving rational functions resulted in high 
efficiency of the method. 
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The outline of the paper is as follows. Section 2 brings a complete solution to the approxi¬ 
mation problem. Some comments on the algorithmic implementation of the method are given 
in Section 3; some technical details of the implementation are presented in Appendix A. In 
Section 4, several examples are given to show the efficiency of the method. In Appendix B, 
some basic information on the Hahn orthogonal polynomials is recalled. 

We end this section by introducing some notation. For y := {yi,y2, ■ ■ ■ ,yd) ^ we 
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denote |y| := yi + ^2 + • • • + Vd, and ||y || := {yf + y^ + ... + yj)^ . 

For n € N and c := (ci, C 2 , C 3 ) € such that |c| < n, we define the following sets (cf. 
Figure 1): 

On ■= {k = (/ci, k 2 ) £ : 0 < |fc| < n}, 


:= {k = {ki,k2) G : ki > ci, /c2 > C2, |fc| < n - C3}, > 

pc — C) \ QC 


( 1 . 1 ) 



c=(0,0,2) 


c= (0,1,2) 


c= (2,1,3) 


Figure 1: Examples of sets (1.1) (n = 11). Points of the set are marked by white discs, while the 
points of the set - by black discs. Obviously, 0„ = U r(). 


Throughout this paper, the symbol H^ denotes the space of all polynomials of two variables, 
of total degree at most n. 

Let T be the standard triangle in 

T := {(xi, X 2 ) : xi, X2 > 0, xi + X2 < 1}. (1.2) 

For n G N, and k := [ki, /C 2 ) G 0^, we denote, 

/ n\ n\ 

Uy A:i!A:2!(n- |fc|)!’ 

The shifted factorial is defined for any a G C by 

(a)o := 1 ; (a)fc := a{a -|- 1) • • • (a -|- A: — 1 ), k > 1. 

The Bernstein polynomial basis in H,^, n G N, is given by (see, e.g., [5], or [6, §17.3]), 

H£(a;) := ^^^xjix2^(l - |a;|)”"l''l, fc := (fei, ^ 2 ) G x:=(xi,X2). (1.3) 
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The (unconstrained) bivariate dual Bernstein basis polynomials [12], 


DU-;cx)eul keOn, 


(1.4) 


are defined so that 

{DI,B^)^ = 5u,u k,leQn. 

Here 6k,i equals 1 li k = I, and 0 otherwise, while the inner product is defined by 

{f,g)a ■= jjWa{x)f{x)g{x)Ax, 

T 

where the weight function Wa is given by 

Wa{x) := ActXi^X2^{l — |®|)“®, OL := (ai, 02 , as), a^ > —1, 


(1.5) 


( 1 . 6 ) 


with Aa := r(|Q:| + 3)/[r(ai + l)r(a 2 + l)r(a 3 + 1)]. 

For n € N and c := (ci, C 2 , C 3 ) G such that |c| < n, define the constrained bivariate 
polynomial space 

:= {p&Ul : P{x) = xrx^^(l - \x\r • Q{x), where Q G nt|e|} • 

It can be easily seen that the constrained set of degree n bivariate Bernstein poly¬ 

nomials forms a basis in this space. We define constrained dual bivariate Bernstein basis 
polynomials, 

Zl^")(.;a)Gn([’2, fcGH);, (1.7) 


so that 


D 


(n,c) 


, Bf) = 6k,I for k, I G 


( 1 . 8 ) 


where the notation of (1.5) is used. For c = (0,0,0), basis (1.7) reduces to the unconstrained 
basis (1.4) in H^. Notice that the solution of the least squares approximation problem in the 
space Hn ’ ' can be given in terms of the polynomials ’ . Namely, we have the following 
result. 

Lemma 1.1 Let F be a function defined on the standard triangle T (cf. {1.2)). The polyno- 

ic 2) 

mini Sn G H),, ’ , which gives the minimum value of the norm 


is given by 


\\F-Sn\\L2 ■= {F-Sn,F-Sn)l, 

Sn=Yl (F,Dt"^) Bt 


(1.9) 




Proof. Obviously, Sn has the following representation in the Bernstein basis of the space 

tt(^>2) . 

rrn 

Jn,cfi 

ken^ 


Sn= 


Bl 
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On the other hand, a classical characterization of the best approximation polynomial Sn is 
that {F — Sn,Q)a = 0 holds for any polynomial Q G (see, e.g. [4, Thm 4.5.22]). In 

particular, for Q = D^’^\ we obtain 

Hence, the formula (1.9) follows. □ 


The coefficients E^{Q.,c,n) in the Bezier form of the dual Bernstein polynomials. 


dI^’^\x;cx) = Y, Et(a,c,n)B^ix), k G H);, 


( 1 . 10 ) 


play important role in the proposed method. Using the duality property (1.8), we obtain the 
following expression for the coefficients of the above expansion: 


E^{a,c,n) = 


{ 


^(n,c) 


,D 


{n,c) 

I 



( 1 . 11 ) 


In a recent paper [11], an efficient algorithm was obtained for evaluating all these coefficients 
for fc, Z G with the computational complexity 0{n^), i.e., proportional to the total number 
of these coefficients. See Section 3.1 for details. 


2 Polynomial approximation of Bezier triangular surfaces with 
constraints 


In this paper, we consider the following approximation problem. 
Problem 2.1 Let R„ be a rational triangular Bezier surface of degree n, 


Rn(®) := 


Y ^kTkBlix) 
Qn(®) _ k£Sn 

~ Y ’ 


£c G r, 


with the control points r^ G and positive weights G M, G ©n- Find a Bezier triangular 
surface of degree m, 

Pm{x) := Y PkBk{x), xeT, 
with the control points pk G satisfying the conditions 


Pk = 9k for k G r^, (2.1) 

gk S being prescribed control points, and c := ( 01 , 02 , 03 ) G being a given parameter 
vector with |c| < m, such that the distance between the surfaces Rn and P^n, 


d{Rn, Pr 


reaches the minimum. 


jjWa{x)\\Rn{x) - Pyn{x)fdx, 

T 


( 2 . 2 ) 
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Remark 2.2 Remember that continuity conditions for any two adjacent triangular Bezier 
patches are given in terms of several rows of the control net "parallel" to the control polygon of 
their common boundary (see, e.g., [6, Section 17]). Therefore, constraints (2.1) are natural, in 
a sense (cf. Fig. 1). In Section 4, we give several examples of practical usage of this approach. 


Clearly, the Bezier triangular patch being the solution of Problem 2.1 can be obtained in 
a componentwise way. Hence it is sufficient to give a method for solving the above problem 
in the case where R„ and are scalar functions, and are numbers. 

All the details of the proposed method are given in the following theorem. 


Theorem 2.3 Given the coefficients and positive weights Wk, k E of the rational 
function 

UkVkB^ix) 

D _Qn(®) _ fceOn 

Rn(^) •— / \ — ^ ; (2.3) 

Y. 


fce©r! 


the coefficients pk of the degree m polynomial 


minimising the error 
with the constraints 
are given by 

where 


Pm{x) := Y PkBkix), 


Rn Pm; Rn Prriion 


Pk = Qk for keVl 


Pk = Y [ 1 ) B^{a,c,m){ui - vi), fc G fl, 

iGoe. 


c 

mi 


ui := 


h€©n 


n\ ( n + m 


1 


VI := 


(|q:| + 3)2m 


hj \ h + I 

E 




JJ(aj + l)hi+ii gh, 


\i=l 


(2.4) 


(2.5) 

( 2 . 6 ) 

(2.7) 


with hs := m — \h\, l^ := m — jZj, and 

rr Bffi^(x) 

JJ \o{x) ^ ^ 

T 

The symbol E^{a,c,m) has the meaning given in (1.10). 
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Proof. Observe that 


l|Rn-Pm|lL = l|W-S^||i^ 

where 

w := - Jm, Tm ■■= •= E 

ker^ 

the notation being that of (1.1). Thus, we want Sm to be the best approximation polynomial 
for the function W in the space Ilm^. Its Bezier coefficients are given by 


where we have used Lemma 1.1. We obtain 


Dr 

(Rn,iir)o= E ^hrni^^Bl 


^ E 

■- E 

h^Bn 


m\ (n + m 


h + l 


-1 


1 


jn+m \ 


dt 


m\ f n + m 


h + l 


-1 


tfi+i, 


where we use the notation (2.8). Further, using equations (1.3) and (1.5), we obtain 

(T„,Bn„= E nkiBt.Br)^ 

r-\h\-\l\ 




Hence, the formula (2.7) follows 


(«! + l)/ii+Zi {.012 + l)/i2+L (“3 + l)2m- 

, ^^\h) V I J dal + 3)2m 


□ 


Remark 2-4 In general, the integrals (2.8) cannot be evaluated exactly. In Section 3.2, we 
show that they can be efficiently computed numerically up to high precision using an extension 
of the method of [9]. 

In the special case where all the weights a;*, i G On, are equal, the rational function (2.3) 
reduces to a polynomial of degree n, so that the problem is actually the constrained polynomial 
degree reduction problem (see, e.g., [18]). Evaluation of the integrals is then a simple task. 


3 Implementation of the method 

In this section, we discuss some computational details of the polynomial approximation of the 
rational Bezier function, described in Section 2 (see Theorem 2.3). 
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3.1 Computing the coefficients E^(a.,c,m) 

We have to compute all the coefficients c,m) with k, I £ It has been shown [11] 

that they can be given in terms of 

with M := m — |c| and p := a + 2c, where are the unconstrained dual 

Bernstein polynomial of total degree M (cf. (1-4)). See Eq. (3.2) for details. Obviously, 
e^(/r, M) = e^(/r, M). The following algoritm is based on the recurrence relations satisfied by 
M), obtained in the paper cited above. 

Algorithm 3.1 (Computing the coefficients E^{a,c,m)) 

Step 1 Let M := m — \c\, p := a + 2c. 


Step 2 For h = 0,1,..., M - 1, 

^2 = 0 , 1 ,... , M — h; 

compute 



(— 1)^1 (|/i,| + 3)m 
M!(/ ri + 2)q 


M-h 

C* hi{l2] IJ.2, fJ-3, M 

i=0 


h), 


(3.1) 


where 0 = (0,0), I = {h,l 2 ), hi{t;a,b, N) are the Hahn polynomials (cf. (B.l)j, 

and 

(/ii + 2 )m . ^ „ 

\ (|p|-W + 2 )M-q’ 

/ (2^ + IpI — /^1 + 1)(/^1 + 2)M-i(|p| + M + 2>)i ^ ^ 

. *!(/^3 + l)i(|/^| ~+ * + l)M-q+l 

next put Cq ■= Cl . 


Step 3 For ki = 0,1,... ,M - 1, 

1° for k 2 = 0,1,... , M — ki — 1, 
h = ki,ki + 1,... ,M, 

^2 = 0,1,... , M — /i, 
compute 

gfc +«2 (^[CTi(fc) - ai(Z)]ef - o- 2 (fc)ef "^2 + o- 2 (/)ef_^J /<Joik), 

where k = {ki,k 2 ), I = {h,h)) '^2 '■= (0,1), and where fort := (ti,t 2 ) we define 


fjo(t) := (|t|-M)(t2+Ai2+l), o-2(t) := t2(|t|-Ai3-Al-l), fTi(t) := f7o(t)+<T2(t), 

next put •= 

2 ° for h = ki + 1, ki + 2,..., M, 
h = 0,1,... , M — li, 
compute 

^k+v^ := ([ri(fc) - Ti(Z)]ef - T2{k) + to(Z) + T2(Z) /ro(fc), 
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where k = (/cijO), I = (^1,^2); Vi := ( 1 , 0 ), and fort := (ii,i2) the coefficients Tj{t) 
are given by 

ro(t) := (|t|-M)(ii+//i + l), T2{t) := ii(|t|-//3-M- 1 ), ri(t) := To{t)+T2{t); 

next put el^^^ := 

Step 4 For k, I ^ compute 

E^{cx, c, m) := t /14 Vi , (3.2) 

where c' := (ci,C2), and 


3 

[/ := (|q;| + 3)2|c| ]^(«i + l)2ci> 
2=1 



As noticed in Remark B.l, the sum in ( 3 . 1 ) can be evaluated efficiently using the Clenshaw’s 
algorithm, at the cost of 0 {M — li) operations. 


3.2 Computing the integrals Ij 

The most computationally expensive part of the proposed method is the evaluation of the 
collection of integrals (2.8). For example, for n + m = 22, if c = (0,0,0), there are 276 
two-dimensional integrals to be computed. It is obvious that using any standard quadrature 
would completely ruin the efficiency of the algorithm. Moreover, if any of the parameters a* 
{i = 1,2,3) in (1.6) is smaller than 0 and the corresponding constrain parameter q equals 
zero, then the integrands in (2.8) are singular functions, and standard quadratures may fail 
to deliver any approximations to the integrals. 

Therefore, for evaluating the complete set of integrals (2.8), we introduce a special scheme 
which is based on the general method [9] for approximating singular integrals. The proposed 
numerical quadrature is of the automatic type, which means that the required number of 
nodes is adaptively selected, depending on the complexity of the rational Bezier function, so 
that the requested accuracy of the approximation is always achieved. Most importantly, the 
algorithm is extremely effective in the considered application. In the example given at the 
beginning of this subsection (n -|- m = 22), the time required to compute the whole collection 
of 276 integrals is only twice^ longer than the time needed to approximate a single separate 
integral of a similar type. 

First, we shall write the integral (2.8) in a different form which is better suited for fast 
numerical evaluation. Observe that bivariate Bernstein polynomials (1.3) can be expressed in 
terms of univariate Bernstein polynomials. Namely, we have 

Bf{x) = Bf^{xi)Bf^~^\x2/il - xi)), j := (ii,j2), x := (xi,X2), 

where Bj^{t) := 0 < i < M, are univariate Bernstein polynomials. Further, 

the bivariate weight function Wa (see (1.6)) can be expressed as 

Wa{x) = Aa Va2+a3,ai (xi) Va3,a2 

^ Based on the Maple implementation of the algorithm. If the collection consists of 990 integrals (n+m = 42), 
the computation time increases by only 50% (compared to the case of 276 integrals). The detailed report from 
the efficiency test can be found at the end of Appendix B. 



where Va,/ 3 {t) := (1 — is the univariate Jacobi weight function. Hence, the integral (2.8) 
can be written as 


"1 (*) 

I-i = I Wa{x) ^ dx2dg;i 



0 JO 


U){x) 


= Ar 


= Ac 


f 


fl 




'*^Q:2+«3 + 1,«1 ('®) '^Ol3,a2{t) 


dt ds 


N 


Va,b{t) / Vc^is) 


1 


00 *{s, t) 


ds dt, 


(3.3) 


.3 / Jo 

where we denoted N ■.= n + m, 

a = a{j) := as + iV - \j\, 
c = c(ji) ■.= a 2 + a 2 , + N 

and 


i=0 j=0 

Note that the computation of values of the integrand is now much more effective, because the 
coefficients Wi of the function w* (1 < i < n) in (3.5) do not depend of the inner integration 
variable s. The main idea is, however, to compute the values of u* only once (at a properly 
selected set of quadrature nodes), and obtain a tool for fast computation of the integrals (3.3) 
for different values of a, b, c, and d, i.e. for different values of j. 

For arbitrary fixed t G [0,1], define the function 


--Oi2+j2, 1 

(3.4) 

= ai+Ji, j 


n—i 

(3.5) 


iptis) ■■= UJ*{s,t) 


-1 


(3.6) 


ft is easy to see that we can write 


Ij = 


with 

where we use the notation 


|J(a,6,ch), 

^>(f) := J(c,(i,V't)> (3-7) 

■■= [ {1 - x)°‘x^f{x)dx. 

Jo 

The functions 'ipt and are analytic in a closed complex region containing the interval [0,1] (it 
is proved in Appendix B). This implies that (cf. [16, Chapter 3]) they can be accurately and 
efficiently approximated by polynomials given in terms of the (shifted) Chebyshev polynomials 
of the first kind. 


Mt 

Ms) ~ SmM) - 1 ), 

i=0 

M 

$(f)«5M(t) := J]'7;7*(2t-1), 

/=o 


0 < s,f < 1, 


(3.8) 
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where M may depend on ji; and the prime denotes a sum with the first term halved. Once 
the above expansions are computed (this can be done in a time proportional to Mt log(Mt) 
and Mlog(M)), the integrals can be easily evaluated using the following algorithm 

that was proved in [14]. 

Algorithm 3.2 (Computing the integral J(a,/3;5), S being a polynomial) 

Given numbers a, /3 > —1, let r := [3 — a, u := a + /I + 1. Let Sm be a polynomial defined by 

M 

Sm{x) = - 1 ). 

i=0 

Compute the sequence dt, 0 < i < A4 + 1, by 


dM+i = dM '■= 0 ) 


di—l ■ — 


2rdi + {i - u)di+i - -g 
i + u 




Output: J{a, /3; Sm) = B ■ (^70 — rdo + udi), where B := r(a + l)r(/3 + l)/r(a + /3 + 2). 


By the repeated use of the above very fast scheme, we may efficiently approximate the 
whole set of integrals Ij for j G The remaining technical details of the adaptive 

implementation of the proposed quadrature and the complete formulation of the integration 
algorithm are presented in Appendix A. 


3.3 Main algorithm 

The method presented in this paper is summarized in the following algorithm. 

Algorithm 3.3 (Polynomial approximation of the rational Bezier triangular surface) 

Given the coefficients r^ and positive weights LOk, k G Qn, of the rational function (2.3), the 
coefficients p^ of the degree m polynomial (2.4), minimising the error (2.5), with the con¬ 
straints (2.6), can be computed in the following way. 

Step 1 Compute the table {E^{oL,c,m)'\f, i^Q^^ by Algorithm 3.1. 

Step 2 Compute the table |A}.,-pQC by Algorithm A.l. 

Step 3 For k G T^, put pk := gu- 
Step 4 For k G 0 ( 7 , compute pk by (2.7). 

Output: Set of the coefficients pk, k ^ 0^- 


4 Examples 

In this section, we present some examples of approximation of rational triangular Bezier 
patches by triangular Bezier patches. No theoretical justihcation is known for the best choice 
of the vector parameter a in the distance functional (2.2) if we use the error function 

A(a;) := ||R„(a;) - Pm(«)|| (4.1) 
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to measure the quality of the approximation. On the base of numerical experiments, we claim 
that cx = ( —usually leads to slightly better results than the ones obtained for 
other "natural" choices of parameters, including the usually preferred cx = (0,0,0) (meaning 
Wa{x) = 1). The computations were performed in 16-decimal-digit arithmetic. In the imple¬ 
mentation of Algorithm A.l, we have assumed e = 5 x 10“^® in (A.3), and used the initial 
values M* = = 32. 


4.1 Example 1 

Let Rg be the degree 6 rational triangular Bezier patch [8, Example 2], 


^6{x) 


uJkrkBlix) 

feG06 

Y ^kBl{x) 

fcG06 


X eT, 


(4.2) 


T being the standard triangle (1.2), and the control points and the associated weights 
being listed in Table 1. We let cx = (—^,—^,—^), c = (1,1,1), and constructed the degree 5 


Table 1: Control points {upper entries) and the associated weights eok {lower entries) of 
the surface (4.2), with k = {ki,k 2 ) € Og. 


ki \ k2 

0 

1 

2 

3 

4 

5 

0 

(6,0,2) 

(5,0,3) 

(4, -0.5, 3.5) 

(3,-0.2, 4) 

(1.5, 0.5, 2) 

(0.4, 0.4,1) 


0.8 

0.3 

1.8 

1.2 

0.8 

0.2 

1 

(5.2,1,3) 

(4.5,1,3) 

(3, 0.6, 4) 

(2,0.9, 3) 

(1.2, 1,2) 

(0.4, 0.8, 0.6) 


1 

0.4 

0.8 

2.4 

1.3 

0.9 

2 

(4.5, 2, 5) 

(4, 2.2, 4) 

(3,2,3) 

(2,1.2, 2) 

(0.8,1.5,1.5) 



0.5 

1 

1 

1.8 

0.8 


3 

(4,3,6) 

(2.5, 2.5, 5) 

(1.5, 2.8, 4) 

(1,2,3) 




0.3 

2 

1 

0.9 



4 

(3.5, 3.5, 4) 

(2.5, 3, 5) 

(1.5, 3.5, 3) 





1.5 

0.6 

1.2 




5 

(3, 4.2, 2) 

(2,4,2) 






0.8 

0.5 





6 

(2,5,1) 

1 






best approximating polynomial patch 






P5( 

X) ■■= Y Pk^Ux), 

X eT, 





fce05 




under the restriction 

Pk = Qk for fc G Tg, where 





:= {k = 

{ki,k2) : ki = 

0 , or k 2 = 

II 

O 

o' 

5}, 


and the set of points gk, k G Tg, is obtained in the following way. As well known, the boundary 
of the patch (4.2) is formed by three degree 6 rational Bezier curves. The least squares degree 
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5 polynomial approximation to each of these rational curves, with the endpoints preservation, 
is constructed using an extension of the method of [14], described in [13] (the input data: 
m = 5, a = /3 = —k = I = 1, notation used being that of [13]). Now, the set of points gu 
is the appropriate collection of all control points of the three resulting Bezier curves. 

We have repeated the computations for cx = (0,0,0) (with a = /3 = 0, in [13]), obtaining 
slightly worse results. The maximum errors maXjcgT ^(®) (cf. (4:.!)) of the obtained results 
(see Fig. 2) are about 50% less than those reported in [ 8 , Table 1]. 






Figure 2: Constrained degree 5 polynomial approximation of the degree 6 rational triangular Bezier 
surface, with c = (1,1,1). Upper part: Rational surface Re{x) and the approximating surface P 5 (a;) 
with a = Lower part: The error A(a;) plots corresponding to a = (0,0,0) and 

CK = (—^, —— 5 )) respectively. The maximum errors are 0.16 and 0.13, respectively. Notice that the 
original surface and the approximating surface agree at the corner points. 


4.2 Example 2 

Let R* be the composite rational surface. 


R*(®) 

where for C £ {R,Y}, 


Rf(z/)) y := {l-\x\,xi-X2), X eTji, 

R^(2), z := {x 2 - xi,l-\ x\), X £Ty, 


R^iw) 


^k^kBliw) 

keBs _ 

^kBl{w) 

fcG05 


w eT, 


(4.3) 


(4.4) 
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T being the standard triangle (1.2), and 


Tr ■. = {x = {xi,X2) ■■ Xi>X2>^, |®| < 1}, 

Ty .={x = {xi,X 2 ) '■ X 2 > xi>Q, |a:| < 1}. 

The control points and the associated weights of the rational patches (4.4) can be 
found at the webpage http://www.ii.uni.wroc.pl/~pwo/programs.htnil. The surface (4.3) 
is shown in Fig. 3 (the left plot). 

Now, we show how to obtain the degree m polynomial approximations of the rational 
subpatches, which form a C^-continuous composite surface. 

1° Let be the triangular Bezier patch of degree m approximating the rational patch 
without constraints, i.e., for c = (0,0,0). Let be the control points of the patch P^. 
2° We approximate the rational patch by the triangular Bezier patch P^ of degree m, 
with constraints of the type c = ( 2 , 0 , 0 ), where the points G T^ are chosen so that the 
C^-continuity is obtained (cf. [ 6 , Section 17]): 

9(0,i) • P(i,0) ’ * 0) 1) ■ ■ ■ ) 

5 (1,0 •= Pji+m + iP{i+i,o) - i = 0 ,1 ,... , m - 1 . 

The results, obtained for m = 5 and m = 6, with a = (—^,—^,—^), are shown in Fig. 3. 
It can be observed that approximation of the rational composite surface (4.3) by two jointed 
polynomial patches of degree m = 5 (the middle plot) resulted in some visible differences. 
Increasing the degree of the approximating polynomials to m = 6 (the right plot) already 
gave a very satisfactory result. 



Figure 3: The composite rational Bezier surface (4.3) {left) and the C^-continuous composite poly¬ 
nomial surfaces of degree (5,5) {middle) and ( 6 , 6 ) {right). 
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Appendix A: The adaptive algorithm for computing the integrals 

h 


We start with proving that the functions (3.6), t G [0,1], and $ (3.7) are analytic in a closed 
complex region containing the interval [0,1]. The assertion is clearly true in the case of = 
uj*{z,t)~^, as the bivariate polynomial uj* has no roots in [0,1] X [0,1]. Similarly, for any 
s G [ 0 , 1 ], the function z e->■ uj*{s, z)~^ is analytic in a rectangular region [—u, 1 + cj] x [—it, o'], 
where o' > 0 does not depend on s. Thus, if s G [0,1], then 


f ui*{s, z) ^dz = 0 

Jc 

for any closed contour C C [—o', 1 + o'] x [—o', o']. Consequently, if q;,/ 3 > —1, then 

(1 — s)“s^u;*(s, z)“^ds^ dz = J {1 — 3 )°" ^ J uj*{s, z)~^ 


ds = 0 . 


IC \ Jo 


Therefore, by Morera’s theorem (see, e.g., [1, Chapter 2.3]), the function ‘h(z) = J{a, j3,'tpz) 
is also analytic in [—o', 1 + o'] x [—o', a]. 

The polynomials SMt Sm in (3.8), which approximates the functions and $, are 
determined to satisfy the interpolation conditions 


SMkisj) = uj*isj,tk) \ 0<j<Mk, 

5 ' m ( 4 ) = J{c,d]SMk), 


0<k<M, 


where, for simplicity, we denote = Mt ^., and the interpolation nodes are given by 


1 1 

Sj = —I— cos-. 

^22 Mk' 


1 1 fevr 

tu = —I— cos —. 

2 2 M 


(A.l) 


In such a case, the coefficients 7 )and 7 ; in (3.8) are given by 



2-5, 


'i,Mk 


Mk 


^cos^. 


Mk 


1=0 


Mk 


11 


2 — 5i^m 
M 


M 


k=0 


COS 


Ikn 

aF’ 


0 < i < Mk, 
0<1<M, 


(A.2) 


where 5j^k is the Kronecker delta, the double prime means that the first and the last term of 
the sum are to be halved. The sets of coefficients (A.2) can be very efficiently computed by 
means of the FFT with only 0(Mfclog(Mfc)) and 0(Mlog(M)) arithmetic operations (cf. [7] 
or [4, Section 5.1]; the authors recall that the FFT is not only fast, but also resistant to 
round-off errors). The presented approach is very convenient from the practical point of view 
because if the accuracy of the approximation (3.8) is not satisfactory, then we may double 
the value of Mk (or M) and reuse the previously computed results. The expansions (3.8) are 
accepted if 


Mk 

E 





maxjl, max 

0<i<3 * 


< 16e 


and 


M 


E It 

i=M-3 


max |l. 


max | 7 i 
0<i<3 


< 256e, 


(A.3) 
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where e is the computation precision. 

Here is the complete algorithm for efficient approximation of the whole set of integrals Ij 
for j G ^n+m- The functions (parameters) a, b, c, and d are defined in (3.4). 

Algorithm A.l (Numerical computation of the set of integrals Ij, j G ^n+m) 

Let M := M*, where M* is an arbitrary integer greater than 7. 

Phase I. For k G {0,1,... , M} do the following Steps 1-6: 

Step 1. Compute t^ according to (A.l), and compute Wi{tk) in (3.5) fori G {0,1,..., n}. 
Step 2. Let := M^, where Mf, is an arbitrary integer greater than 7. 

Step 3. Compute the values uj*{sj,ti^)~^ for j G {0,1,... , M^}, where Sj is given by 
(A.l). 

Step 4- Using the FFT, compute the coefficients (0 <i < M^) defined in (A.2). 
Step 5. If the first condition of (A.3) is not satisfied, then set := 2Mk, compute 
the additional values uj*{sj,tk)~^ for j G {1,3, 5,... , — 1}, and go to Step f- 

Step 6. Compute the set of quantities W[tk,ji] := J {c{ji),d{ji by applying 

Algorithm 3.2, for ji G (ci, ci + 1,..., A — C 2 — C 3 }, where N = n + m. 

Phase II. For ji G {ci, ci + 1,..., A — C 3 — C 2 } perform the following Steps 7-9: 

Step 7. Compute the coefficients 7 ; (0 < I < M) defined in (A.2), by means of the 
FFT, using the stored values W[tk,ji], 0 < k < M, in place of j(^c{ji),d{ji)', Sm^) ■ 
Step 8. If the second condition of (A.3) is not satisfied, then set M := 2M, and repeat 
Steps 1-6 for A; G (1, 3,5,... , M — 1}. 

Step 9. For j 2 G {c 2 , C 2 + 1,... , A — C 3 — ji}, compute the integrals 

^3 = kh,32) '■= T(a(j), 6 (i 2 );^M) 

using Algorithm 3.2. 

Output: Set of the integrals Ij for j G 


Remark A.2 In Steps 4 and 7 of the above algorithm the coefficients (0 < z < M^) or 7 ; 
(0 < Z < M) are recalculated each time the value of Mk or M is doubled. Such a procedure 
is advised if we use a system (like, e.g.. Maple or Matlab) equipped with a fast built-in FFT 
subroutine. If we are to program the FFT summation algorithm by ourselves, it should rather 
be done in such a way that practically all results computed for a previous value of or M 
are reused (cf., e.g., [7]). 

In Table 2 we present the results of the efficiency test, where the proposed quadrature (imple¬ 
mented in Maple) is compared to the Maple built-in integration subroutine. We have used the 
Bezier surface form Example 4.1 (n = 6 ), and set the parameters m and c to several different 
values, to obtain collections of integrals of different sizes (equal to |fl()_|_^|). The experiment 
was performed in the 64-bit version of Maple 16 on the computer equipped with the 3.7GHz 
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Table 2: Comparison of the computation times of the Maple library function and the proposed 
adaptive quadrature (Algorithm A.l) in the case of several collections of integrals (2.8). The number 
of integrals which are to be computed equals 


\^n+m\ 

computation time (in seconds) 

Maple library function 

the proposed method 

1 

0.064 

0.30 

3 

0.19 

0.30 

10 

0.64 

0.32 

28 

1.75 

0.37 

91 

6.34 

0.43 

276 

22.9 

0.59 

990 

FAILURE 

0.89 


17 processor. All parameters a* in (1.6) were set to 0 (the efficiency of the proposed method 
does not depend on a, but the Maple built-in integration subroutine works most efficiently 
with this selection). 

We have to keep in mind that Maple is an interpretative programming language with 
a pretty slow code interpreter. Therefore, the 4.7 times longer computation time of our 
quadrature, compared to the computation time of the Maple library function, in the case 
of 1-element collection of integrals is in fact an excellent result. The last collection of 990 
integrals {n + m = 42) was too difficult to be computed by the Maple built-in subroutine (in 
14-decimal digit arithmetic, assumed during this test). 


Appendix B: Hahn orthogonal polynomials 

The notation 

f ) ■ ■ ■ j Or I \ (ai)fc • • • (Qr)fc k 

^ • ^^kKh)k---{bs)k 

is used for the generalized hypergeometric series (see, e.g., [2, §2.1]); here r, s € z, ai,bj G 
C, and (c)k is the shifted factorial. The Hahn polynomials (see, e.g., [10, §1-5]) 

hi{t) = hi{t- a, b, M) := (a + 3 F 2 ^ ^ 0 ’ 

where / = 0,1,..., M, a, b > —1, and M € N, satisfy the recurrence relation 

hi+i{t) = Ai{t,M)hi{t) + Bi{M)hi-i{t), I > 0] hoit) = 1] h-i{t) = 0, (B.2) 

with the coefficients 

Ai{t,M) := Ci{2l + s — l) 2 t — Di — El, Bi{M) := —Di Ei^i, (B-3) 

where s := a -|- 6 -|- 1, Ci := ( 2 / -|- s -|- 1 )/[(Z -|- s)( 2 Z -|- s — 1 )], Di := Ci l{l + M + s){l -|- b), and 
El := (/ + a + l)(M-/). 

Remark B.l A linear combination of Hahn polynomials, SN{t) := 'Yld=Qlihi{t',a,b,M), can 
be summed using the following Clenshaw’s algorithm (see, e.g., [4, Thm 3.2.11]). Compute the 
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sequence Vb, Fi,..., K +2 from Vi := 7 * + Ai{t; M)Vi+i + Bi+i{M)Vi+ 2 , i = N, N - 1,... ,0, 
with V^+i = yAf +2 = 0, where the coefficients and Bi{M) are defined by (B.3). 

Then S]si{t) =Vq. 
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